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ABSTRACT 

fH | We present a detailed investigation of the effects that a direct interaction be- 

tween Dark Energy and Cold Dark Matter particles imprints on the Halo Mass 
Function of groups and clusters of galaxies. Making use of the public halo catalogs of 
the CoDECS simulations, we derive the Halo Mass Function for several different types 
of coupled Dark Energy scenarios both based on the Fricnds-of-Friends algorithm and 
on the Spherical Overdensity halo identification for different values of the overdensity 
threshold A c . We compare the computed Halo Mass Functions for coupled Dark 
Energy cosmologies with ACDM as well as with the predictions of the standard 
analytic fitting functions. Our results show that the standard fitting functions still 
reproduce reasonably well both the Friends-of-Friends and the Spherical Overdensity 
Halo Mass Functions of interacting Dark Energy cosmologies at intermediate masses 
and at low redshifts, once rescaled to the characteristic amplitude of linear density 
perturbations of each specific model as given by as- However, we also find that 
such apparent degeneracy with as is broken both by the high-mass tail and by the 
redshift evolution of our Halo Mass Functions, with deviations beyond ~ 10% for 
most of the models under investigation. Furthermore, the discrepancy with respect 
to the predictions of standard fitting functions rescaled with the characteristic value 
of as shows - for some models - a strong dependence on the spherical overdensity 
threshold A c used for the halo identification. We find that such effect is due to a 
significant increase of halo concentration at low redshifts in these models, that is 
however absent in the majority of the cosmological scenarios considered in this work. 
We can therefore conclude that the universality of the Halo Mass Function is violated 
by cosmological models that feature a direct interaction between Dark Energy and 
Cold Dark Matter. 
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1 INTRODUCTION 

Despite its longstanding success in explaining a wide range 
of astrophysical and cosmological observations, the standard 
cosmological model based on a cosmological co nstant A as 
the d r iver of the observed accelerated expansion (iRiess et al.l 
Il998l : iPerlmutter et al. | GUI; ISchmidt et al.l 1 19981 ) and on 
Cold Dark Matter (CDM) particles as the source of struc- 
ture formation processes, is far from being accepted as a 
satisfactory description of the Universe. The lack of a firm 
physical interpretation of the two main dark constituents 
of the cosmic energy budget, Dark Energy (DE) and CDM, 
represents in fact an open issue that is presently motivat- 
ing an enormous effort in the community both concerning 



the development of alternative theoretical models and the 
implementation of observational tests capable to falsify the 
standard scenario. 

Furthermore, the severe fine-tuning probl ems that 
chara cterize the cosmological constant (see e.g. IWeinberel 
1989) and the increasing number of astrophysical observa- 
tions suggesting a possible ten sion with the predictions of 
the ACDM paradigm (s ee e.g. iBov lan-Kolchin et al. 201 ll: 
Lee fc Komatsul l20ld : iThomas et~ai1 1201 ll : lLovell et all 
20111 ) have stimulated the development of a wide vari- 
ety of alternative scenarios. These range from simple dy- 
namical DE models base d on a minimal l y coupled scalar 
field as e.g. Quintessence l|Wetterichlll988l ; iRatra fc Peebles! 



2 W. Cui, M. Baldi, S. Borgani 



1988) or k-essence (lArmendariz-Picon et al. 2001 ) to in - 
teracting DE models jWetterichl 1 19951 : lAmendolal l200d ), 
to modifications of General Relati vity at cosmological 
scales (see e. g. Hu fc Sawicki 120071 ). Warm Dark Mat- 



ter scenarios jBode et all 2001 ), or Clustering DE models 
ijCreminelli et alJl20ld )~ All these models should therefore 
be tested by a direct comparison of their characteristic fea- 
tures with observations. This requires to develop accurate 
predictions of the effects of all these different scenarios on 
various observable quantities, from the background evolu- 
tion down to the highly nonlinear regime of cosmic structure 
formation 

In the next decade, a large number of observational 
ventures, such as the ESA approved EUCLID^ mission 
l|Laureiis et al.l 1201 ll ). will enormously improve the con- 
straints on the minimal set of parameters that character- 
ize the standard ACDM cosmological scenario. However, the 
quest for a physical interpretation of the dark components 
of the Universe demands for an extension of these observa- 
tional probes beyond such minimal parameters set, allowing 
for a richer phenomenology of the dark sector with respect 
to what envisaged by the standard model. In fact, although 
clearly necessary, the continuous improvement of the obser- 
vational constraints on the standard cosmological parame- 
ters might not be sufficient to detect possible failures of the 
ACDM scenario, thereby preventing a deeper understand- 
ing of the physical properties of the dark universe. This is 
mainly due to the fact that the observable deviations from 
the expectations of the standard model are in many cases 
totally or partially degenerate with one or more of the stan- 
dard cosmological parameters, which makes it particularly 
difficult to disentangle the two effects. 

It is therefore a crucial task for a full exploitation of 
present and future observational probes to investigate to 
which extent the characteristic features of a wide range of 
alternative cosmological scenarios might be hidden beyond 
a simple shift in one or more of the standard ACDM cos- 
mological parameters. In particular, in the present paper we 
will study in great detail the halo mass function (hereafter 
HMF) of groups and clusters of galaxies in the context of 
a significant number of interacting Dark Energy scenarios 
and we will investigate their apparent degeneracy with the 
amplitude of linear density perturbations encoded by the 
power spectrum normalization <T8. Indeed, a precise calibra- 
tion of the HMF represents the basic ingredient to exploit 
the pote ntial of galaxy clu sters as tracers of cosmic growth 
(see e.g. jAllen et alj|201ll . for a recent review). 

Within the standard ACDM paradigm, a number of 
HMF studies have been carried out through the years, 
based on N-body simulations covering progressively larger 
dynamic ranges, with the purpos e of calibrating fittin 



functions for a universal HMF (e.g. ISheth fc Tormenl 
Jenkins et all l200ll ; iReed et al.l 120031 ; IWarren et al] 
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2006; 



Lukic et alj|2007f l. or to characterize subtle deviations from 
this universality (e .g. IReed et al.l 120071 ; iTinker et all |200S| ; 
ICrocce et al"]|2010h . Extensions of these studies to include 
the effect of baryons on the HMF have been made possible 
only recently by the advent of cosmological hydrodynami- 
cal simulations covering an adequate range of halo masses 



jRudd et al.ll2008l ; IStanek et al.ll2009l ; ICui et al.ll201lh . It is 

however clear that the quest of exploring cosmological mod- 
els beyond the standard one, requires enlarging the range 
of models for which the HMF calibration needs to be car- 
ried out, and whose universality with respect to variations 
of the parameters relevant for each sp ecific class of mod- 
els ne eds to be assessed. For instance, iBhattacharva et al.l 
Hon]) analyzed the HMF for an extended suite of simula- 
tions of quintessence models, and confirmed violation of uni- 
versality at the ~ 10 per cent level for the range of masses 
and redshift covered by their simulations. 

A semi-analytic study of the HM F for coupled DE mod - 
els has been recently carried out by iTarrant et al.l (|20 1 if ) . 
Using the expressions for a universal ACDM HMF, they 
computed the variations in the cluster counts induced by 
the modification in the spherical collapse model induced by 
the DE coupling to CDM. 

In this paper, we make use of the publi cly available cat- 
alogs of the CoDEC^I N-body simulations (|Baldill2011cl '), to 
analyse in detail the effect that different coupled DE models 
have on the evolution of the HMF. We will directly compare 
the number counts of massive halos identified both through 
a Friend-of-Friend (FoF) algorithm and through a Spherical 
Overdensity (SO) algorithm at different overdensities with 
the standard fitting functions calculated assuming different 
values of as- Former analysis of the HMF for N -body simu- 
lation s of C oupled DE mo dels were presented in lBaldi et al.l 
l|20ld ) and lBaldil j2011dl ). Their analysis, which was based 
on simulations covering a narrower dynamic range and pa- 
rameter space than those analysed here, and only considered 
halos identified from the FoF algorithm, showed that HMF 
in Coupled DE models can be described well by the ACDM 
HMF after suitably rescaling erg. This implied that the signa- 
ture of coupling in the observational analysis of the cluster 
mass function could be completely masked by a modified 
power spectrum normalization. In our analysis we will ver- 
ify whether this degeneracy between DE coupling and power 
spectrum normalization can be broken at the level of accu- 
racy allowed by these simulations, by following the redshift 
evolution of the HMF. 

The paper is organized as follows. In Section 2 we de- 
scribe the Coupled DE models considered in our analysis. 
Section 3 is devoted to the description of the corresponding 
simulations and to the halo identification procedure. In Sec- 
tion 4 we will present the results of our analysis of the HMF 
and their interpretation in terms of variation of halo con- 
centration induced by the presence of the DE coupling. We 
will summarize our main results and draw our conclusions 
in Section 5. 



2 THE COUPLED DARK ENERGY MODELS 

Interacting DE models have been widely discussed 
and inves tigated in the l iterat u re in the last years 
see e.g. lAmendolal |2000|. 120041 ; iFarrar fc Rosenl 1 2007 ; 



Pettorino fc Baccigalupil 120081; lAmendola et alJ 12008 
Kovama et al.l 120091 ; iHonorez et alj l2010l ; iBaldil l2011d 



1 |http777www. euclid-ec.org/ 
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-1.0 
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-OL(f) 
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V{<j>) = Ae 


-OL(f) 
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V(4>) = Ae 
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0.08 


0.15 





-0.992 


0.967 


EXP008e3 


V(4>) = Ae 


-OL(f> 


0.08 


0.4 
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-0.982 


0.895 


SUGRA003 


V(4>) = A<t> 


-Q e <* 2 /2 


2.15 


-0.15 





-0.901 


0.806 



Table 1. The list of cosmological models considered in the CoDECS project and their specific parameters. 



iTarrant etall l201ll ; IClemson et~ai] l20Tlf ) and a new de- 
tailed presentation of these cosmological scenarios would 
be superfluous in the present paper. We therefore refer 
the interested reader to the above mentioned literature for 
a thorough description of coupled DE (cDE) models and 
for the derivation of their main equations. We limit our 
discussion here to the definition of the notation and of the 
conventions adopted for the specific cDE models considered 
in our analysis. 

In the present paper, we will consider the set of cDE 
models pres ently include d in the CoDECS suite of N-body 
simulations iBaldii l201ld ) - the largest set of cosmological 
simulations to date for interacting DE cosmo l ogies - that 
have been p resented and discussed in IBaldii l|2011bl ) and 
IBaldii (|2011ch . These are fiat cosmological models where the 
role of DE is played by a dynamical scalar field <f> with a self- 
interaction potential V(4>) exchanging energy-momentum 
with the CDM fluid through an interaction term denned 
by the following set of equations: 



p c + 3Hp c 

p b + 3Hp b 
p r + 4Hp r 

SH 2 



(Pr + pc + pb + P4>) 



(1) 

(2) 

(3) 
(4) 

(5) 



where the subscripts b, c, and r indicate the baryonic, CDM, 
and radiation components of the universe, respectively. In 
Eqs. (11 loll an overdot represents a derivative with respect 
to the cosmic time t, H = a/ a is the Hubble function, and 
Mpi = I/VStvG is the reduced Planck mass. The source 
terms at the right hand side of Eqs. (1 1 12l! represent the 
interaction between DE and CDM, where the dimension- 
less coupling function /3 c (</>) sets the coupling strength while 
the sign of the quantity /3 c ((j>)(j> defines the direction of the 
energy-momentum flow between the two components. The 
energy exchange determines a time variation of the CDM 
particle mass, according to the equation: 



d In M c 



(6) 



dt V 3 

which can be derived from Eq. ([2]). 

In the present work, we will consider two possible 
choices for the coupling function /3 C (4>), defined as: 

p c (J>) = (3 e M , (7) 

namely a constant coupling = 0) and an exponen- 



tially gro wing coupli n g (/3i > 0). The latter case, first pro- 
po sed by lAmendolal l|2004j ) and subsequently investigated 
by IBaldii l |2011dh . allows for lar ger values of the present 
coupling strength /3q as compared to constant coupling 
models, since the impact of the interaction on the back- 
ground expansion history of the universe and on the Cos- 
mic Microwave Background anisotropics is strongly sup- 
pressed by the time evolution of the scalar field <j)- Fur- 
thermore, we will consider two distinct choices also for the 
scalar self-inte raction potential V{4>\ namely an exponen- 
tial potential llLucchin fc Matarresei Il985l : IWetterichlll988l : 
iFerreira fc Jovcelll998l ~ 



VO) = Ae' ' 4 ' (8) 
and a SUGRA potential ijBrax fc Martinll 1999ft : 

V(cj>) = 4f°e* 2/2 , (9) 

where for simplicity the field <j> has been expressed in 
units of the reduced Planck mass in Eqs. (|8l9fl . The main 
phenomenological difference between these two potential 
functions resides in the existence of a global minimum at a 
finite <f> value for the SUGRA potential, while the exponen- 
tial potential is monotonically decreasing to zero for <f> — > 00. 
The presence of a global minimum in the SUGRA potential 
allows for an inversion of the scalar field motion and for a 
consequent change of sign - in case of a constant cou pling 
/3 C - of the quantity /3c4> , as discussed in lBaldil (|2011bj| ) (see 
also ITarrant et all l201ll ). Due to such inversion, the DE 
equation of state parameter w$ = p^/ ' p$ shows a "bounce" 
on the cosmological constant "barrier" w$ = — 1, for which 
this class of models ha s been dubbed the "Bouncing cDE 
scenario" (Baldi 2011b). For the specific model considered 
in the present work, the "bounce" happens at relatively 
recent epochs, Zi nv ~ 6.8, and has significant consequences 
on the ev olution of lin ear and nonlinear perturbations 
(see again IBaldii l2011bh . The effect of the coupling on 
the background evolution of the universe is to allow for a 
phase of Early Dark Energy which goes under the name 
of -MDE (0-Matter Dominated Epoch, see lAmendolal 
|2000|) or G-0-MD E (Growing-0-Matter Dominated Epoch, 
see iBaldii l2011dh for models with constant and variable 
couplings, respectively. Such scaling behavior of the DE 
density determines a different expansion history of cDE 
models with respect to a ACDM cosmology with the same 
cosmological parameters, which represents one of the most 
characteristic features of cDE scenarios and that is correctly 
taken into account in the numerical implementation of the 
CoDECS simulations described in the next Section. All 
the features and the parameters of the different models 
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investigated in the present work are summarized in Table[T] 



Parameter 



Value 



The effect of the coupling on the growth of linear density 
perturbations is described by the evolution equations for 
linear fl uctuations 5 b , c for baryons and CDM, respectiv ely 
(see e.g. lAmendola|[2004 iPettorino fc BaccigalupilhoOST ): 



8 C = -2H 



Hy/6 

S b = -2H5 b + 4ttG [p b S b + p c 5, 



8 C + AtvG [p b S b + p c S c T c ] , (10) 
(11) 



where for simplicity the field dependence of the coupling 
function I3 c {<f>) has been omitted. The factor r c is defined 
as T c = 1 + 4/3;? (</>) /3 and represents a fifth-force acting on 
CDM perturbations, while the term 2/} c (</>)0/v6 is an extra- 
friction arising as a consequence of momentum conservation. 
The combination of these two effects determines the devi- 
ation of the evolution of linear density perturbations from 
the standard ACDM case. At the nonlinear level, the same 
two extra terms ap pear in the accele ration equation of cou- 
pled particles fsee lBaldi et al.ll2010l . for a derivation of the 
acceleration equation) : 



v c 



Y> GM C {4>)V C + ^ GM b 



n 



(12) 



where r c ^ b is the physical distance of the target coupled par- 
ticle from the other CDM and baryonic particles, respec- 
tively. However, in this case the relative orientation of the 
particle's velocity v c and of the local gradient of the grav- 
itational potential V<f> introduces an additional degree of 
freedom and plays an important role in the evolution o f non- 
linear structures in the context of cDE scenarios (see iBaldil 
l2011al . for a detailed discussion of nonlinear effects in cDE 
models) . 

Clearly, for any chosen coupling /3 c ((j>) and potential 
slope a in a cDE model, it is always possible to reconstruct 
an effective equation of state w e a(z) for a minimally coupled 
scalar field cosmology that provides the same expansion his- 
tory. In this respect, the specific background evolution of 
any cDE model does not represent a "smoking gun" capa- 
ble to uniquely identify the presence of an interaction in the 
dark sector, as the same expansion history might be pro- 
vided by a suitably tuned standard Quintessence or phantom 
DE fluid. On the other hand, any pair of cDE models with 
different coupling functions but with identical cosmological 
parameters at z — (except for the value of the equation of 
state) will necessarily have different background expansion 
histories due to their different </!>-MDE scaling solutions in 
matter domination (since the energy fraction of the Early 
Dark Energy component during </>-MDE is p roportional to 
/3 2 ((j>), see e.g. lAmendol^l200d ; lBaldll2011dl ). Therefore, in 
order to investigate solely the impact of the DE-CDM cou- 
pling on the growth of structures one would need to compare 
any given cDE model with a specific standard Quintessence 
or phantom DE scenario having the same expansion history. 
However, this would necessarily imply to choose a different 
reference model for every different coupling function j3(<f>). 
In our study, instead, we are interested in comparing a range 
of different cDE models with one single minimally coupled 
reference scenario, given by the concordance ACDM cosmol- 
ogy, all having the same WMAP7 parameters at z = 0, by 
including all the effects (i.e. modified background expan- 



H 70.3 km s" 1 Mpc" 1 

fi CDM 0.226 

tt DE 0.729 

o-g 0.809 

Q b 0.0451 

n s 0.966 



Table 2. Parameters of the reference ACDM model which ia 
used to generate initial conditions for all the cDE simulations 
(see text). 



sion history, enhanced growth of linear density perturba- 
tions, particle mass variation, and scalar fifth-force). 



3 THE SIMULATIONS AND HALO 
IDENTIFICATION 

3.1 The N-body simulations 

For our analysis we will make use of the public data of the 
CoDECS simulations fealdil l2011<J ) - the largest suite of cos- 
mological simulations to date for cDE cosmologies - which 
include all the specific models listed in Table [T] The pub- 
lic data of the CoDECS project have already been used to 
investigate the characteristi c features of cDE s cenarios in 
terms of z-space distortions (IMarulli et al.ll2011f). of the ex- 
pected Weak Lensing signatures (IBevnon et alJl201lT) and o f 
the abundance of "Bullet-like" systems |Lee fe Baldill201lh . 

In the present work, we will consider in particular the 
L-CoDECS runs, that follow the evolution of 1024 3 CDM 
and 1024 3 baryon particles in a periodic cosmological box 
of 1 comoving /i _1 Gpc aside. The sim ulations hav e been 
carried out with the modified version bv lBaldi et all (|2010h 
of the widely u sed parallel TreePM N-body code GADGET 
l|Springel|[2 005). which implements all the specific features of 
cDE models described in Sec. [21 The L-CoDECS simulations 
have a mass resolution at z = of m c = 5.84 x 10 10 /i _1 Mq 
for CDM and m b = 1.17 x 10 10 /i _1 M Q for baryons. 

In order to cover a large dynamic range over which to 
study the effect ou coupled DE with an affordable computa- 
tional cost, our simulations account for the different forces 
acting on the uncoupled DM and coupled baryonic matter 
components, without including the hydrodynamic descrip- 
tion of baryons. As such, our simulations are designed to fol- 
low the evolution of two populations of collisonless particles, 
which f eel different g r avitat i onal forces. W e refer t o the anal- 
yses by iRudd et al.l (|200Sh : IStanek et ail (|2009h ; ICui et all 
(2011) for discussions on the effect of hydrodynamics on the 
HMF. 

Initial conditions for all the simulations have been 
generated b y perturbing a homogeneous "g lass" particle 
distribution (|Whitelll994l ; iBaugh et al"] Il995h according to 
Zel'dovich approximation 1 Zel'dovichl [1970) in order to 
produce a random-phase realization of the same initial 
matter power spectrum, corresponding to th e one com- 
puted by the public Boltzmann code CAMB (|Lewis et al.l 
|2000T ) for a ACDM cosmology with parameters consistent 
with the "WMAP7 o nly Maximum Likelihood" results of 
iKomatsu et al.l (|201lh . which are listed in Table [2] 
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Since the same random seed has been adopted for the 
initial conditions generation of all the L-CoDECS runs, all the 
different cosmological models share exactly the same parti- 
cle displacements at zcmb ~ 1100. The initial conditions for 
each specific simulation are realized by scaling such displace- 
ments to the starting redshift of the numerical integration 
Zi = 99 with the specific growth function obtained by nu- 
merically solving Eqs. (| 101 11 p for each individual cosmology. 



3.2 The Halo Identification 

The two most common methods for halo identification in 
N-body simula tions are based on the Friend-of-Friend (FoF) 
algorithm (e.g. iDavis et al.lll985 | ) and on the sph erical over- 
density (SO) algorithm (e.g. iLacev fc Colelll994h . The FoF 
halo finder has only one parameter, b, which defines the 
linking length A as A = b ■ I, where / = n -1 ^ 3 is the mean 
inter-particle separation, with n the mean particle number 
density. In the SO algorithm there is also only one free pa- 
rameter, namely the mean density A c p cr it(z) contained 
within the sphere over which the halo mass is computed, 
with pcrit(z) being the critical cosmic density at redshift z. 
Note that our definition of the spherical overdensity A c is 
given with respect to the critical density p cr it instead of the 

which means that our 



background matter density p„ 
A c is Q m times the commonly used overdensity parame- 
ter A. Each of the two halo finders has its own advantages 
and shortcomings (see more d etails in IJenkins etall l2001: 
IWhitdl2003 ; lTinker et al.ll200Sl . and references therein). The 
difference of halo mass and HMF resulting from the two 
methods ha v e been discussed in several analysis (see e.g . 



White! [200: 



More et 



oo5 



. .Reed et al l l2003t 120071 : ICohn fc White] 120081; 
20111 ). 



In this paper, we applied both methods to identify halos 
and to compare the respective HMF for the CoDECS set of 
N-body simulations. First, the standard FoF algorithm has 
been applied, with a linking parameter b = 0.2, only over 
the CDM particles in the simulations as primary tracers of 
the matter distribution. Subsequently, each baryonic parti- 
cle has been attached to the FoF group of its closest CDM 
neighbor. The gravitationally bound substructures of each 
FoF halo have also been identified b y means of the SUB- 
FIND algorithm (see more details in ISpringel et"aH l200ll ; 
iDolag et al.ll2009h . The position of the most bound particle 
of each substructure defines then the center of the subhalo. 

As a second step, the SO halo catalogues are built using 
a fast SO algorithm (see also lCui et al.ll201lf ). In this SO al- 
gorithm, the centers of spherical regions are identified with 
the positions of the most bound particles within the main 
subhalo of each FoF group, as identified by SUBFIND. Mean- 
while, the region encompassed by R200 (i.e. the spherical re- 
gion enclosing a mean overdensity 200 times larger than the 
critical density) around each FoF halo is also checked, and 
all the substructures identified by SUBFIND, that lie outside 
this region, are considered as new independent FoF groups. 
The most bound particle of each such new independent ha- 
los is then used again as the center of a new SO region. 
For each independent halo identified with this method we 
progressively increase the radius of a sphere centered on the 
most bound particle until a specified mean internal density 
contrast A c p cr it(z) is reached. The mass Ma c within this 



spherical region of radius _Ra c is then 
4 , 

Ma c = -TvR Al: A.cPcrit(z) . (13) 

Since each halo is firstly identified starting from a FoF 
algorithm, it inherits some FoF disadvantages. A well known 
potential problem with FoF is that there are situations in 
which two halos are connected through a bridge of particles. 
Our check of all the substructures lying outside -R200 should 
reduce this ef fect, especia l ly for the most massive halos. As 
discussed by iReed et al.1 (|2007l ). this effect becomes more 
important at high redshift and for poorly resolved low-mass 
halos. Since our analysis mainly concentrates on very mas- 
sive halos, which we select to be resolved by at least 200 
particles, we expect the bias induced by using FoF parent 
groups to be relatively small. In order to verify this on a few 
test cases, we also used the SUBFIND results to check the 
mass ratio of all the subhalos within each FoF group, with 
respect to the main halo. Whenever a prominent subhalo is 
found whose mass is at least 0.7 times the mass of the main 
halo, we treat it as a distinct halo. We verified that such ha- 
los are extremely rare after we remove subhalos lying outside 
-R200, at least over the mass range M > 10 13 /i _1 Mq consid- 
ered in our analysis, so our results are not influenced by this 
effect. Finally, since the groups identified by the FoF algo- 
rithm have by definition no overlapping, we do not include 
in our identification of SO halos any res triction to prevent 
such overlapping (see lTinker et all 120081 for a discussion on 
halo overlapping). 

W e do not apply the c ommonly used FoF mass correc- 
tion o f IWarren et all i|200rj ) for our FoF halos in this paper. 
This correction is meant to regulate the halo mass, which 
is overestimated by the FoF algorithm when halos are sam- 
pled with too few particles. This correction is expected to 
be small for our selected halo sample, due to the minimum 
threshold of 200 particles per halo that we adopted. 



4 RESULTS 



Under simplifying assumptions, iPress fe Schechterl (|l974j) 
were the first to provide a theoretical description for the 
abundance of collapsed halos as a function of their mass. 
Thanks to the calibration on N-body simulations cov- 
ering wide dynamic ranges, the HMF for the standard 
ACDM cosmology can now be accurately described by an- 
alytical models and fit ti ng functions, a s prov i ded by e.g. 
Sheth fc Tormenl (I1999T); IJenkins et all d200ll); IReed et al.l 



(|2003h : IWarren et all (|2006l ); iTinker et al.1 l|2008T ). In most 
applications, the HMF is expressed in terms of the so-called 
multiplicity function by: 



M 



dn(M, z) 



p (z) diner" X (M, ; 



(14) 



where n(M, z) is the number density of halos with mass M 
at redshift 2 and po{z) is the background density at redshift 
z. In Eq. (|14|) . a 2 (M,z) is the variance of the linear den- 
sity field, extrapolated to the redshift z at which halos are 
identified, which is given by 



\M,z) = 



2tt 2 



k 2 P(k)W 2 (k,M)dk. (15) 
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Figure 1. The FoF halo mass function (HMF) for the different cDE simulations at z = 0, 1 and 2 (left, central and right panels, 
respectively). In all cases b = 0.2 is assumed for the linking parameter used in the FoF halo identification. In the upper part of each panel 
th e symbols show the re sults from our simulations, while curves are the corresponding predictions from the FoF HMF of Eq. ll7l calibrated 
by Ijenkins et al ] 4200 lF l for a ACDM cosmology and rescaled to the value of erg appropriate for each cDE model. In all panels, black 
crosses and black curves are for the results of our ACDM reference simulation. The lower part of each panel shows the difference with 
respect to ACDM results. Here symbols show the difference between cDE and ACDM simulations, while curves with different line-styles 
show the difference between the corresponding fitting functions. 



In Eq. (I15|) . D(z) is the growth factor of linear density per- 
turbations, P(k) is the linear power spectrum of density per- 
turbations and W(k, M) is the Fourier-space representation 
of a real-space top-hat filter of radius R, which on average 
encloses a mass M = 4ttR 3 po(z)/3, given by: 



W{k,R) 



(kRY 



[sm(kR) - kRcos(kR)} . 



(16) 



Using a series of N-body simulations under the standard 
ACDM cosmology and ap plying the Fo F algo rithm with 
linking parameter b = 0.2, Ijenkins et al.l l|200ll ) (J01 here- 
after) provided a simple fitting formula as a function of the 
variance a, which is independent of redshift: 



/(a) = 0.315 exp(-| lna -1 + 0.61| 3 8 ) 



(17) 



This fitting function is valid over the range —1.2 ^ 
lncr -1 ^ 1.05, which corresponds to halo masses 10.3 ^ 
log(M//i _1 M Q ) < 15.6 in our ACDM simulation. We use 
this formula as our benchmark reference fit for the FoF 
HMF. We adopt instead the fitting function of lTinker et al.l 
(|2008h (T08 hereafter) for the SO HMF. This fitting func- 
tion f(a) is tuned by the four parameters A, a, b, c according 
to: 



m = A 



+ 1 



(18) 



These four parameters are functions of two independent vari- 
ables, namely the redshift z and the overdensity A (for more 
details, see Eqs. 5-7 and B1-B4 in T08). As already men- 
tioned, T08 adopted an overdensity A computed with re- 
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Figure 2. The ratio between the FoF halo counts in different mass intervals from the cDE models and from the reference ACDM 
model from simulations, compared to the same quantity computed from the corresponding fitting functions from Eg. 1211 (see also text). 
Deviations of this quantity from unity indicates that simply rescaling trg in the fitting function of Eq |17l does not correctly account for 
the difference between cDE and ACDM simulation results. The correspondence of simulated models with symbols and line-types is the 
same as in Fig. [T] 



spect to the background matter density. In order to use this 
fitting function with our definition of the overdensity A c 
(i.e. with respect to the critical density of the Universe), we 
simply change A to A c in the functional dependence of the 
four parameters (A,a,b, c) just multiplying it by Q m (z). 

Throughout this paper, we adopt the differential HMF 
dn/d log to make comparisons between theoretical predic- 
tions and simulation results. With a simple transformation 
of Eq. (|14p , the theoretical expression for dn/d log M is given 
by: 



dn 



dlogA/ 



f{a)po ding 1 
M dlogM ' 



(19) 



where f(cr) is given by Eqs. (|17|l and Q5) . 

In order to compute the differential HMF dn/d log M in 
a simulation of volume L 3 , one needs to measure the number 
of halos AN in a given logarithm mass bin A log M : 



dn 



M AN 



dlogM L 3 A log M ' 



(20) 



where we assigned the characteristic mass M for a given 
mass bin as the mean mass computed over all halos belong- 
ing to that bin. For our ha lo sample, we used narrow mass 
bins with A log M = 0.2. iLukic et all i|2007l ) had studied 
the error introduced by binning data, and showed that this 
error is negligible with such a narrow bin. However, the lim- 
ited statistics of high-mass halos within the simulation box 
makes the determination of the HMF quite noisy in the high 
mass end, especially when marrow mass bins are used. To 
reduce this effect, we merge mass bins containing less than 
20 objects into the adjacent lower mass bin. Each mass bin 
is then weighted proportionally to the number of clusters 
it contains. Due to this specific treatment of the mass bins 
at large masses, the bin width can be different for different 
simulations. Therefore, when we compare the number of ob- 
jects in such last bins, we rescale the cluster counts within 
each of them to the bin width in units of A log M = 0.2. 



4.1 The FoF HMF in coupled dark energy models 

We show the FoF HMF of the different cDE models included 
in the CoDECS project derived from the L-CoDECS simula- 
tions in the upper panels of Fig. [1] Here, the symbols show 
the results from simulations, while the different line-styles 
represent the HMF predicted by the J01 fitting function 
given by Eq. (|17[1 for ACDM cosmologies with the same as 
value as the different cDE realizations (see Table [TJ. In the 
lower panels of each plot of Fig. [T] we show with symbols the 
ratio of the halo number density of the various cDE simula- 
tions over the corresponding results for the standard ACDM 
simulation. In the same panels, the curves show the same ra- 
tio, but obtained using the model predictions from the HMF 
fitting functions, each computed for the appropriate value of 
us- In fact, since the HMF is described only by the variance 
of the linear density field cr 2 (M, z), the HMF fitting func- 
tions of the ACDM cosmology are often assumed to predict 
also the HMF of other non-standard models by a simple 
renormalization of the linear perturbations amplitude a$. 

At a first sight, by looking at Fig. Q] such assumption 
holds also for cDE models considered in the present paper, 
with the predicted HMFs for the as values computed for 
each model through linear perturbations theory fitting the 
cDE simulations results reasonably well, at least at z — 
(left panels). This is true for both the standard cDE mod- 
els with constant coupling, EXP001-3 (upper plots), and for 
the other two cDE models with different potentials or cou- 
pling functions, EXP008e3 and SUGRA003 (lower plots). 
Quite interestingly, the accuracy of such rescaling degrades 
at increasing redshift, with deviations appearing at z = 1 
(central panels), which becomes more apparent at z = 2 
(right panels). 

In order to better appreciate the difference between re- 
sults from cDE simulations and fitting functions, we intro- 
duce a parameter T(M, z), which is expressed as 



T(M, z) = 



NB2(M,z)/NZgh M (M,z) 
N f D %{M,z)/N^ DM {M,z) 



(21) 



where N(M, z) is the number density of halos within the 
mass bin M at redshift z, and the lower index of indi- 
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cates the model, while the upper index indicates either the 
theoretical {fit) or the simulation (sim) results. With this 
definition, the quantity T(M, z) describes the difference be- 
tween simulations and fitting functions, after rescaling each 
to the corresponding ACDM prediction. In this way, we ac- 
count for sampling effects, due to the limited halo statistics 
especially in the high-end of the simulated HMF, and for 
possible inaccuracies in the fitting functions. Accordingly, a 
significant deviation from unity of the T(M, z) parameter 
indicates a lack of precision in the rescaling of the HMF 
fitting function for the cDE models. 

In Figure [5] we show the mass dependence of T(M,z) 
at the same redshifts considered in Fig. [T] for the different 
cDE models. At z = 0, T is very close to unity at small halo 
masses, thus implying that the J01 fitting function can be 
used to capture the simulation results quite accurately. The 
only exception is represented by the "Bouncing cDE model" 
SUGRA003 (red line), that shows a ~ 5 per cent deviation 
at the low-mass end of our HMF, with a decreasing trend at 
higher masses. Deviations from unity of the T parameter for 
this model change their sign at higher redshift, where they 
become progressively larger. At z — 2 simulation results for 
the SUGRA003 model show a deviation from ACDM simu- 
lation results that is ~ 30 per cent smaller than predicted by 
the corresponding J01 fitting function, for the highest sam- 
pled halos masses. An opposite trend is instead found for 
the variable-coupling model EXP008e3, with T ~ 1.4 for 
the highest masses sampled at z = 2. As for the standard 
models with constant coupling (EXP001-3), the value of T 
shows smaller deviations from unity, which also in this case 
increases with mass and redshift. As expected, the smallest 
deviations are found for the models with the smallest value 
of the p coupling (EXP001). 

In general, Figure [2] clearly shows that, even for the 
FoF HMF, the apparent degeneracy between the coupling 
and as, which holds to first approximation in the low-mass 
end at low redshift, is broken in the high-mass end by an 
amount which increases with redshift, with the strength of 
the coupling, and whose sign also depends on the shape 
of self-interaction potential for the dynamical scalar field 
which determines the coupling. 

4.2 The SO HMF in coupled dark energy models 

The FoF algorithm does not assume any geometry for the 
halo, unlike the enforced sphericity of SO. For this reason, 
the mass function for halos identified with a SO algorithms 
can be more directly compared with the observed mass func- 
tion of galaxy clusters and groups, whose mass is generally 
measured, or inferred from mass proxies, within some radius 
encompassing a fixed overdensity A c with respect to the cos- 
mic critical density. In this section, we investigate the SO 
HMF at three reference overdensities A c = 200, 500, 1500. 

In Figure Owe show the same SO HMF properties at 
A c = 200 as shown for the FoF HMF in Figure [T] above. 
From the upper part of the panels of Figure [3] we see that 
at a first sight the SO HMF for different cDE models can 
also be described by a simple renormalization of as in the 
T08 fitting function of Eq. [THJ The ratios to ACDM in the 
lower panels of Figure[3]are also qualitatively similar to what 
already shown for the FoF HMF in Figure □ All the cDE 
models based on an exponential potential present more mas- 



sive halos than ACDM at all the three redshifts considered, 
while the "Bouncing cDE model" SUGRA003 is consistent 
with ACDM at z = 0, but shows an excess of massive ha- 
los at higher redshifts, in line with the results discussed by 
iBaldil l|2011bl ). Also in this case, we show the T parameter 
as a function of halo mass in Figure U Differently from what 
shown for the FoF halos in Figure [2] the T parameter of the 
three standard cDE models (EXP001-3) are quite close to 
unity with a decline < 10 per cent at all redshifts and at all 
masses (except at the very high mass end). The T parame- 
ter for the EXP008e3 model shows an increase with redshift, 
which is in any case still within ~ 10 per cent. On the other 
hand, the T parameter for the SUGRA003 model is above 
unity, by ~ 20-40 per cent, at z — 0. Its value is larger than 
what shown for the FoF case of Fig. [5] while dropping below 
unity at z = 2 by ~ 10-20 per cent. In general, these results 
suggest that the as rescaling applied to the T08 fitting func- 
tion provides a better fit to the SO HMF at A c = 200 for the 
standard cDE models than the same rescaling does when ap- 
plied to the J01 fitting expression for the FoF HMF, except 
for the bouncing cDE model SUGRA003. 

The SO HMFs at A c = 500 are shown in Figure [5] 
At this higher overdensity, the SO HMFs from simulations 
are still fit by the theoretical predictions as in Figure [3] with 
the ratio displayed in the lower panels also showing a similar 
quality of the fit as for the previous cases. From Figure [S] 
we see that the T parameter for the standard cDE models 
has larger deviations from unity than for the A c = 200 case, 
especially at higher redshift. As for the SUGRA003 model, 
this parameter is confirmed to have a decreasing trend with 
redshift, starting with a value at z — which is larger than 
the corresponding one found for A c = 200. 

In general, all these trends are amplified as we con- 
sider progressively larger values of A c . In fact, at our high- 
est overdensity A c = 1500, the fitting function by T08 does 
not provide a good fit to the simulation results anymore, 
even for the ACDM simulation as shown in Figure [7] Nev- 
ertheless, the ratios to corresponding ACDM results have 
a similar variation with halo mass as in the previous fig- 
ures. Since the T parameter is the ratio between simulation 
and fitting-function results, rescaled to the corresponding 
ACDM predictions, it should allow to remove this offset. 
Figure[S]indicates that this is indeed the case, with the value 
of T still being within ~ 20 per cent from unity at z = 
for all the exponential potential models, while at z — 2 we 
find deviations from unity of the order of ~ 30 per cent for 
the most extreme EXP003 model. On the contrary, for the 
SUGRA003 model T increases above unity by more than 60 
per cent at z = 0, and drops to ~ 1 at z = 2. Some devia- 
tion from a universal shape of the HMF has been detected 
also for uncoupled DE scenarios based on a classical scalar 
fiel d with no dire c t inte ractions to matter, as discussed by 
e.g. lCourtin et all (|201lT ) . However, the higher amplitude of 
the deviation and its peculiar mass dependence - with the ef- 
fect being progressively more pronounced for halos of higher 
masses - make the imprint of cDE models on the HMF quite 
peculiar and potentially distinguishable from the one of an 
uncoupled model. 

As a general result from the analysis of the SO HMF, 
we point out that rescaling the mass function in cDE models 
with the as value produces results that, at least for A c = 
200, are more accurate than for the FoF HMF. However, 
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Figure 3. The same as in Figurc[T] but for the HMF com puted with the sphe rical overdensity (SO) algorithm at A c = 200. In this case, 
curves are obtained from the fitting function of Eg. I18l bv iTinker et al.l (j2008T l - computed for the appropriate overdensity. 
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Figure 4. The same as in Figure[2] but for the HMF com puted with the sphe rical overdensity (SO) algorithm at A c = 200. In this case, 
curves are obtained from the fitting function of Eg. U8l by |Tinker et alj l|200St ). computed for the appropriate overdensity. 
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Figure 5. The same as in Figurc[T] but for the HMF com puted with the sphe rical overdensity (SO) algorithm at A c = 500. In this case, 
curves are obtained from the fitting function of Eg. I18l bv iTinker et al.l {2003), computed for the appropriate overdensity. 




13.5 14.0 14.5 15.0 13.5 14.0 14.5 13.2 13.4 13.6 13.8 14.0 

log(M/h"'M e ) logCM/VX) log(M/h"'M ) 



Figure 6. The same as in Figure[2] but for the HMF com puted with the sphe rical overdensity (SO) algorithm at A c = 500. In this case, 
curves are obtained from the fitting function of Eg, U8l bv lTinker et alj l|200St ). computed for the appropriate overdensity. 
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this rescaling does not provide accurate predictions of the 
cDE HMF at higher overdensity and larger redshift. These 
results highlight that degeneracy between a non-vanishing 
coupling and erg can in principle be broken by observational 
measurements of the cluster mass function over a sufficiently 
large mass and redshift range. 



4.3 Halo Concentration 

As shown in the previous section, deviations of the SO HMF 
for cDE models from the ACDM one increase at higher over- 
density. This result suggests that this difference in the HMF 
is due to a change in the timing of halo collapse induced by 
the presence of coupling which, in turn, induces a change of 
the halo density profiles. 

Halos in dissipationless N-body simulations of a ACDM 
model have spherically aver aged density profiles th at are well 
described by the profiles of iNavarro et al l (|l997h . hereafter 
NFW. For an NFW halo of a given mass, the halo den- 
sity profile can be specified entirely by one parameter, the 
concentration, although the relation between concentration 
and mass is characterized by a significant intrinsic scatter 
(e.g. iMaccjo et alll200sl; iGao et al.ll2008l ; iDuffv et al.ll200Sl ; 
IZhaoet alJl2009l : IPrada ,et al.ll201lf ). 

In the following , we use a simple method (see Eqs. 6-9 in 
ISpringel et al . 2008) to compute concentrations for the halos 
identified in our L-CoDECS simulations. According to this 
method, the characteristic NFW overdensity 5 C is a function 
of halo concentration c, and can be expressed in terms of the 
maximum circular velocity of the halo, Vmax , and the radius 
T m ax at which this velocity is attained: 

_ 200 c 3 / V max y 

3 m(l + c) - c/(l + c) ' \H r max ) ' 

(22) 

In Fig. |5J we show the mean concentration c calcu- 
lated from Eq. (|22p . as a function of the FoF halo mass. All 
the cDE models appear to be consistent with the ACDM 
results, with the only exception being represented by the 
SUGRA003 model. For this model the average concentra- 
tion is a factor of ~ 2 larger than in ACDM at z = 0, 
while gradually approaching the ACDM relation at high red- 
shift. This enhanced halo concentration in the SUGRA003 
model could provide a way to account for the apparent over- 
concentration of massive clusters recently found from lens- 
ing analyses (e.g. lOguri et al.ll2009l ; [Umetsu et al.ll201ll ). On 
the other hand, if such a high concentration would hold 
also for typical CDM halos of luminous spiral galaxies (a 
few x lO 12 /!^ 1 A/q), which are not resolved by the simula- 
tions analysed here, this specific realization of the "Bounc- 
ing cDE" scenario could run into strong tension with the 
observed dyna mical properties of spira l galaxies, as recently 
investigated bv lBaldi fc Salucca (|201ll ). The enhancement of 
halo concentration at low redshifts for the "Bouncing cDE" 
scenario also explains why the HMF of this model shows 
deviations from ACDM only at high overdensity thresholds, 
A c = 500 , 1500, while the FoF HMF and the low overdensity 
SO HMF (i.e. for A c = 200) is much closer to the expecta- 
tion from the fiducial ACDM model. 



5 DISCUSSION AND CONCLUSION 

We have investigated in large detail how the halo mass func- 
tion is affected by different possible types of interaction be- 
tween dark energy and dark matter. By me ans of the L - 
CoDECS suite of large N-body simulations i|Baldil [20lTbh 
we have computed the abundance of halos as a function of 
redshift and mass in six different cosmological models (see 
Table 1): the fiducial standard ACDM scenario, three mod- 
els of interacting DE with an exponential self-interaction 
potential and a constant coupling function (EXP001-003), 
one model with an exponentially growing coupling strength 
(EXP008e3), and one specific realization of the Bouncing 
cDE scenario based on a SUGRA self-interaction potential 
(SUGRA003). In doing so, we have assumed as a single ref- 
erence scenario the standard concordance ACDM cosmol- 
ogy with cosmological parameters consistent with the latest 
WMAP7 results, and we have constructed all the other cDE 
models such that they have the same WMAP7 parameters 
at z — 0. This strategy allows us to perform a self-consistent 
comparison of the footprints of cDE models with respect to 
the standard ACDM cosmology. On the other hand, in order 
to distinguish the peculiar signatures of a DE-CDM inter- 
action on the growth of structures, from the effects of any 
other minimally coupled dynamical DE model (such as a 
Quintessence or a phantom scalar field) one would need to 
compare any given cDE scenario to a corresponding dynam- 
ical DE cosmology tuned to have the same expansion his- 
tory. For each of our scenarios, we have then derived both 
the Friends-of-Friends (FoF) halo mass function and a se- 
ries of SO mass functions for three values of the overdensity 
threshold A c , namely A c = 200,500,1500, at the present 
time and at z — 1 and 2. 

Our findings show that the FoF mass function is signif- 
icantly affected by the interaction and displays large devi- 
ations from the expected ACDM halo abundance at z — 0, 
due to the faster growth of density perturbations, with the 
only exception of the Bouncing cDE scenario that shows 
very little differences from ACDM at the present time. This 
is expected due to the peculiar dynamics of the Bouncing 
cDE model that provides the same value of erg as ACDM 
at z — 0. In this respect, our results fully confirm previous 
findings. On the other hand, all the remaining cDE mod- 
els are characterized by larger values of as as compared to 
ACDM, and our investigation shows that the effect of the 
interaction on the mass function is in general highly degen- 
erate with (78, such that the nonlinear halo mass function 
extracted from the simulations can be quite accurately re- 
produced by the standard fitting functions computed for the 
actual value of as attained by each model. Such degeneracy 
is however broken by the different redshift evolution of the 
halo mass function, especially at large masses, in the vari- 
ous cDE models as compared to ACDM, such that the same 
value of as does not accurately reproduce the abundance 
of massive halos at different epochs. This is true also for 
the Bouncing cDE model, which for low overdensity thresh- 
olds appears practically indistinguishable from ACDM at 
z — 0, but shows a significant excess in the expected num- 
ber of massive halos at higher redshifts. On the other hand, 
at higher overdensity thresholds the Bouncing cDE model 
can be easily distinguished from ACDM even at z — 0, as 
the degeneracy with as is broken and the halo mass func- 
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Figure 9. The mean halo concentration as a function of FoF halo mass. The different colors, symbols and line-styles are the same as in 
the previous figures. The shaded region shows ltr scatter in the concentration computed within each mass bin. For reasons of clarity it 
is only shown for the reference ACDM simulation. 



tion shows a significant enhancement over the expected halo 
abundance at large masses. We have investigated the ori- 
gin of this dependence of the deviation from ACDM for the 
Bouncing cDE scenario, and found that the effect is due 
to a significant increase of the average halo concentration 
at low redshifts for such cosmology: starting from the same 
normalization of the concentration-mass relation as ACDM 
at z = 2, the Bouncing cDE model produces concentrations 
that are twice as large as i n ACDM at z = 0. 

In a recent analysis iTarrant et al.l (|201ll ) computed 
the linearly-extrapolated critical overdensity for spherical 
collapse, <5„, for coupled DE models, with the purpose of 
including the effect of coupling in the halo mass function. 
Clearly, including the coupling effect in the computation 
of 8, affects the p redicted halo abu n dance when using 
the expressions by Press fc S chcchter; l|l974f ) (PS) or by 
ISheth fc Tormenl (ll999Tl (ST ) . On the other ha nd, the fitting 



Jenkins et ail |2001) (J01) and iTinker et at 



functi ons by 

(2008) (T08) are not affected by <5* since their universality 
makes them to be only functions of the r.m.s. fluctuation 
amplitude computed at the halo mass scale. While it is 
beyond the scope of this paper to test the predictions of PS 
and ST mass functions by accounting for DE coupling on 
spherical collapse, our results clearly demonstrate that the 
universal fitting functions by J01 and T08 can not be used 
to predict mass and redshift dependence of halo counts 
with an accuracy of 10 per cent or better. This is even more 
true when dealing with non-standard cDE models, such as 
the EXP008e3 and SUGRA003 ones. 

To conclude, we have investigated the effects of interact- 
ing Dark Energy models on the abundance of massive halos 
as a function of redshift and mass, using different methods 
to identify halos and to compute their mass. Our results 
show that a clear degeneracy exists between the coupling 
and the standard parameter as, but that both the redshift 



evolution and the detailed shape of the high-mass tail of the 
halo mass function allow to break such degeneracy at some 
level. Furthermore, we have shown that halo concentrations 
are very mildly affected by the Dark Energy interactions 
when the amplitude of density perturbations is normalized 
at high redshifts (contrarily to what is found for a normal- 
ization at z = 0) with the only exception of the "Bouncing" 
coupled Dark Energy scenario that shows a very rapid in- 
crease of halo concentrations with respect to the fiducial 
ACDM model at recent epochs. Our study therefore pro- 
vides a direct way to test interacting Dark Energy models 
with present and future data on the abundance of massive 
clusters as a function of redshift. A full exploitation of future 
cluster surveys to constrain coupled DE models will however 
require accurately calibrating corrections to expressions of 
the halo mass function, whose universality has been tested 
with simulations only within the ACDM framework. Clearly, 
the observational measurement of the evolution of the clus- 
ter mass function alone could not provide in itself an incon- 
trovertible test for the presence of cDE. This test should be 
complemented by other observational probes of large-scale 
structures. In this framework, the calibration of the HMF 
for cDE models presented here add an important piece of 
information that, in combination with other observational 
cobstraints, will enable to test and, possibly, to falsify the 
cDE scenario. 
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